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Until now, most of our knowledge about the universality class of crystal plasticity has come 
from simulations using discrete dislocation dynamics. These are force-controlled, typically at zero 
temperature, and deal with the creation and annihilation of dislocations phenomenologically. In 
this work, we go beyond these limitations by using phase field crystal simulations in two dimensions 
at finite temperature to extract the avalanche statistics of a simulated crystal under constant shear 
velocity. In addition to the avalanche size and energy distributions we extract the avalanche duration 
distributions and power spectra. All exponents and scaling functions extracted here for the statics 
and dynamics of crystal plasticity, belong to the mean field elastic depinning universality class, 
confirming earlier findings based on discrete dislocation dynamics. 

PACS numbers: 62.20.F-, 61.72.Bb, 89.75.Da, 64.60.De 


Crystal plasticity is in some sense a solid state analogue 
of fluid turbulence, with deformation at micro scales be¬ 
ing both intermittent and spatially inhomogeneous Hi- 
Ill]. These phenomena have been captured realistically 
by a number of approaches, including discrete dislocation 
dynamics (DDD) models [U ITTI - [T7| . continuum models 
[l2|, phase field models [HI [19] and phase field crystal 
(PFC) models |20j . The most prominent open question 
that remains unanswered is the one that unifies the accu¬ 
mulated literature and solidifies the hard-earned knowl¬ 
edge: What is the universality class of crystal plasticity? 

A wealth of critical exponent and scaling function 
information from DDD simulations (in two m Ha 
and three dimensions m) and experiments on slowly 
compressed nanocrystals and microcrystals 0 uni m] 
strongly suggest that crystal plasticity belongs to the 
mean field depinning universality class. Nevertheless the 
issue is by no means settled; for example recent DDD 
simulations in 2D obtained non-mean field results [22] . 

Our main aim in this paper is to uncover the univer¬ 
sal behavior of deforming crystalline matter as it em¬ 
anates from microscopic origins and percolates through 
all scales. To this end we study crystalline plasticity 
with a phase field crystal model [501 155H55] . The ele¬ 
mentary entity in our simulations is a phase field repre¬ 
senting local atomic density, which is appropriately con¬ 
strained to behave as an atomic crystal, and indeed can 
be related to density functional theory |5T|. The phase 
field crystal exhibits elastic, reversible deformation at 
small external loads. It deforms plastically at exceed¬ 
ingly large external loads. Large deformations imprint 
permanent, irreversible change in the lattice structure. 
At sufficient shearing the periodicity is broken and topo¬ 
logical defects emerge in the system. Dislocations travel 
and interact with each other through the lattice forming 
intricate structures. One can observe plastic deformation 
being mediated through intermittent dislocation motion, 
i.e. through discrete slip-avalanches. Here we extract 


several avalanche measures and show that they are dis¬ 
tributed according to power-laws over several orders of 
magnitude revealing long range spatial and temporal cor¬ 
relations. The set of critical exponents we calculate (the 
duration distribution power law exponent for the first 
time) fully supports the mean field depinning picture for 
crystal plasticity. This is in strong agreement with earlier 
2D DDD simulations m- It also agrees with 3 dimen¬ 
sional simulations of dislocation dynamics m, with an¬ 
alytics [iniiis] and experiments PHHIIII]. Our work ad¬ 
dresses the fundamental question of the universality class 
of crystal plasticity and is relevant to the deformation of 
nano-crystals [211 EH EH] and micro-crystals [I10IIQ] as 
the need for miniaturization of devices expands both in 
breath and depth. 

The Phase Field Crystal Model, Sheared: The phase 
field crystal (PFC) model [23] EH 132] describes how the 
local density of atoms changes with time while main¬ 
taining the symmetries and periodicity of the lattice. In 
addition, the elastic interactions of the atoms are also 
captured by the PFC allowing for the elasticity of the 
crystal to be expressed. These characteristic properties 
of the phase field crystal model are distinctly different 
from the phase field model. In a typical phase field model 
the phase field describes the dynamics of interfaces that 
separate dissimilar regions without keeping track of the 
microscopic information inside those regions. In [181119] 
Koslowski et al. developed a phase field model to simu¬ 
late dislocations as interfaces (separating crystal regions 
with different accumulated slip). In studying plasticity, 
however, it is important to capture the microscopic de¬ 
tails such as the dislocations which disrupt the periodic¬ 
ity of the perfect lattice. At the same time, it is impor¬ 
tant to capture the macroscopic behavior as well such as 
the collective motion of the dislocation ensemble. The 
phase field crystal model is particularly successful in do¬ 
ing that in an elegant way [20]. 

The total free energy in the phase field crystal (PFC) 
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TABLE I: Table of exponents. Our results from 2D PFC simulations are shown in the fourth column while our results from 
2D DDD simulations are shown in the third column. Mean field interface depinning values are in the fifth column. Results 
from a full 3D DDD simulation are indicated with an asterisk (*). Results from a 2D DDD simulation with creation and 
annihilation in the steady state are indicated with a plus sign ("*"). Symbol definitions: D{x) stands for the distribution of x, 
Xmax is the maximum of the distribution of x, PS stands for power spectrum, (x) stands for average of x, S is the size of a slip 
avalanche, tavai is its duration, E its energy. V{t) = X/i-i collective dislocation speed, r stands for shear stress, 

Tc is its critical value. Small greek letters are used for critical exponents throughout. 
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is a functional of p{x, t) the local density of the phase field 
(at point X = {x,y) in space and time t). The first term 
in Eq.Q penalizes departures of p{x,t) from periodic¬ 
ity, thus describing a crystal structure as a density wave. 
The last two terms impose a double well potential (to 
lowest order) similar to the Landau ansatz. The reduced 
temperature r is given by {T — Tc)/Tc and controls the 
phase behavior. The material is liquid for temperatures 
higher than a critical temperature while it crystallizes 
for temperatures below T^. Thus, for r > 0 one finds a 
liquid, constant p, phase (because the potential is single 
well) while for r < 0 the stable state is a triangular lattice 
or a striped phase (due to the double well potential). 

The PFC density p{x, t) evolves according to the mas¬ 
sive phase field crystal equation HOI 133]: 
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where a controls the range and j3 the time scale of 
phonon excitations of the crystal |33j . Thermal fluctu¬ 
ations are represented by the stochastic noise p which 
is assumed to be Gaussian with second moment given 
by the fluctuation-dissipation theorem {p{x^ t)p{x', t')) = 
—eW'^5{x — x')5{t — t'). The noise amplitude sets the scale 


of temperature e ^ ksT. This free energy is governed by 
conservative, relaxational, diffusive dynamics that can be 
derived from density-functional theory [32) . 

The PFC solid is a perfect triangular lattice in equi¬ 
librium, but its excitations are phonons and topological 
defects such as dislocations. When the phase field crys¬ 
tal is sheared, it will respond by generating dislocations, 
just as a real crystal. The ability to create and annihilate 
dislocations naturally and easily is one of the advantages 
of the phase field crystal method, compared to DDD. 

We applied a shear strain rate along the x direction at 
the y = 0,Ly boundaries by adding the convection term 
v{y)dp/dx, to the evolution equation Eq. (H- The bound¬ 
ary shear velocity profile, v{y) = (yo = 0 

for -b, yo = Ly for —) is designed to be mainly controlled 
by the velocity at the boundary vg since its penetration 
length X Ly does not affect the results strongly. The 
simulations take place in a square box of sides L^^Ly in 
the x^y direction. The boundary conditions are periodic 
in X and fixed at j/ = 0, Ly i.e. we design the simula¬ 
tion cell such that the crystal wraps around in a; = 0, 
and terminates at ?/ = 0, Ly (without wrapping around). 
That way we can easily apply a fixed shear rate at the 
y = 0, Ly boundaries and allow the dislocations to flow 
unbounded through the x = Q,Lx boundaries effectively 
simulating a larger thermodynamic system than the mere 
dimensions of our basic simulation cell. 

The PFC model has the added value over DDD simula¬ 
tions that it incorporates nonlinear elasticity [34] as well 
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as dislocation creation and annihilation seamlessly with¬ 
out requiring additional phenomenological rules to model 
these number-changing operations. The PFC methodol¬ 
ogy is thus uniquely capable of addressing such questions 
as how strain heterogeneity drives dislocation number 
fluctuations, which in turn couple to plasticity avalanches 
[HH] . Since the PFC model is in essence an atomic simula¬ 
tion, the highly nonlinear interaction at small dislocation 
distances is captured naturally through the crystal lattice 
that mediates it. The same holds true for the particular 
effects of creation and annihilation of dislocations. 

The PFC model handles applied shear velocity na¬ 
tively. DDD models incorporate applied external stress 
naturally. Thus the PFC is suitable to investigate the slip 
avalanches above the critical point (flow stress) in the de- 
pinned state while the DDD below the critical point, in 
the pinned state. In that sense they perform complemen- 
tarily to each other. 

PFC simulations at Finite Shearing Rate and Temper¬ 
ature: We study crystal plasticity as it proceeds intermit¬ 
tently through slip avalanches using the sheared phase 
field crystal (PFC) model [20]. We obtain the main scal¬ 
ing behavior of the distribution of a variety of avalanche 
measures and for several different temperatures (e) and 
shearing rates (vq). We find remarkable agreement be¬ 
tween simulations and analytical mean field theory pre¬ 
dictions of exponents [lOlllllSSj. Our results strongly 
support the critical point picture of plasticity, and sug¬ 
gest new experiments. 

At every time step we obtain the phase field density 
p in a 2D square simulation cell (i.e. = Ly = L) 

through Eq.([^. Large values of p indicates PFC ‘atoms” 
while low field signifies interatomic space, cooperatively 
arranged into a tight crystal (triangular in 2D). By apply¬ 
ing shear along the fixed boundaries y = 0,L the trans¬ 
lational symmetry breaks and dislocations are created in 
an attempt to relieve the high stress accumulated near 
those boundaries. Once the dislocations are created, they 
interact with each other to form pairs and more complex 
structures such as low-angle grain boundaries; individual 
dislocations, dislocation pairs and grain boundaries can 
be seen in the snapshot of the PFC simulations in Supple¬ 
mentary Material. Dislocations may also glide through¬ 
out the entire crystal, allowing for slips to self-organize 
into slip avalanches. We quantify the avalanche activity 
by extracting the total speed of the dislocations, 

N{t) 

Vit) = ^ IF.I, (3) 

i=l 

where N{t) is the number of dislocations in the sys¬ 
tem at time t and Vi is the speed of the i-th disloca¬ 
tion. This measure is similar to the acoustic emission 
signal in Weiss et al.’s single crystal ice experiments (e.g. 
[Tj). Other variants of the avalanche activity measure 
in the literature can track the avalanches as well. For 


example the collective dislocation velocity is defined as 
the dislocation’s burgers vector) 
and represents the strain rate of the crystal [ini[l2ll30|. 
Note that all dislocations, single dislocations or disclina- 
tion pairs that move at high speeds and grain boundaries 
that move slowly, participate equally in the calculation 
of V(t) (see figure in Supplementary Material). 

A perfect unsheared triangular crystal is devoid of dis¬ 
locations and therefore every atom has = 6 nearest 
neighbors (identified with Delaunay triangulation). Con¬ 
versely there exists one dislocation for every atom with 
Hi = 5 or Hi = 7 neighbors since vacancies are not al¬ 
lowed. 

We capture the speed of the dislocations, V (t), through 
the speed of these defect atoms with m ^ 6 neighbors, 
V(t) = ^^detBct(t) |-.|^ where Afdefect(t) is the number of 
defect atoms and Ui is the velocity of defect atom i. 

In order to partition the signal into individual slip 
avalanches we apply a threshold, Fthr, to it for each tem¬ 
perature and shearing rate we simulate. The beginning of 
the avalanche is signified at an instant when the collective 
dislocation speed V{t) intersects upward the threshold 
while its end is when V{t) crosses the threshold down¬ 
ward immediately after. We extract the probability dis¬ 
tribution of the avalanche duration 


^aval ^finish ^starti 


(4) 


size (also called activity fluctuations in the flowing state) 
S = V(t)dt, and energy E = V‘^{t)dt, where 

fstart and tfinish are the starting and ending time of the 
event respectively. In order to see the fluctuations that 
correspond to slip avalanches we applied a threshold 
equal to the average of the signal, V (t ), in each realiza¬ 
tion of total time ftotai, Vthreshoid = F(t)dt, 

We also calculate the power spectrum of V(t): 


PS(uj) 
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(5) 


The power spectrum reveals the frequency content of the 
time series of the collective dislocation speed and it needs 
no thresholding. It is equivalent to extracting the time¬ 
time correlations of the collective dislocation behavior 
and near a critical point it is expected to have no char¬ 
acteristic scale, i.e. be scale-free [571155] . 

For each shearing rate and temperature, we run 48 
different realizations each with different seed for the ran¬ 
dom number generator for the noise p to obtain sufficient 
statistics. This results in tens of thousands of avalanche 
events for each shearing rate vq and temperature e. In 
Fig. [T] we show the event size, duration and energy dis¬ 
tributions while in Fig. [^the power spectra and average 
size versus duration for different shearing rates at the 
same temperature. We find that the distributions follow 
a power law for small event sizes and cut off at larger 
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sizes, with the maximum avalanche size not exhibiting 
a strong dependence on shear rate over the simulated 
range of the parameter vg- We suspect that the reason 
that shear rate does not affect the distribution much is 
that for our systems the system size sets the cutoff of 
the avalanche size distribution. For much larger systems 
the mean field theory predicts that an increase in shear 
rate will reduce the cutoff of the avalanche size distribu¬ 
tion. This can only be seen in systems that are so large 
that the system size is much larger than the correlation 
length of the avalanches given by the finite shear rate. 
In previous work with PFC HO] a different threshold was 
applied to the signal so as to extract avalanches. It re¬ 
sulted in a rate-dependent avalanche distribution cutoff. 
Although both thresholds reveal the same power-law in 
the avalanche distributions we believe the threshold used 
here (equal to the average of the dislocation activity for 
each run) is a more natural way to quantify the fluctua¬ 
tions around the mean dislocation activity. 



FIG. 1: (color online) Probability distributions Dt(tavai) 

(left), Ds{S) (middle) and De{E) (right) of avalanche du¬ 
ration tavai, size S and energy E respectively for different 
shearing rates vo at the same temperature e = 1.6. The prob¬ 
ability distribution of the slip avalanche sizes follows a power 
law with exponent k ~ 1.5, of durations with Kt ~ 2 and of 
energies with ke ~ 1.3. These results are in agreement with 
MFT (Table m. 

Scaling Behavior of the avalanches: The distributions 
of the avalanche size, duration and energy, the power 
spectra and average size versus duration shown at same 
shearing rate, vq = 0.765, and different temperature pa¬ 
rameter values e are shown in figures in the Supplemen¬ 
tary Material. In Figs, [l] and we presented the dis¬ 
tributions of the avalanche size, duration, energy, power 
spectra and average size versus duration at the same tem¬ 
perature parameter values e = 1.6 and different shearing 
rates, vg. Each curve is characterized by a power law for 
several decades and a cutoff at large values (smaller val¬ 
ues for the power spectra) which does not change with 
shear rate. Mean field theory predicts the dependence 
on the shear rate, that should be visible in larger sim- 



FIG. 2: (color online) The power spectrum PS{u}) of 

the collective dislocation speed (main) and the average slip 
avalanche size (S) versus duration tavai (inset) for different 
shearing rates vg at the same temperature e = 1.6. The 
power spectrum scales with the inverse square frequency giv¬ 
ing Xjavz ~ 2. The average size scales with the square of 
the duration, l/avz « 2, for avalanches that are sufficiently 
small not to touch the sample boundaries. (The boundaries 
of the power-law scaling regime of the power spectrum PS{u)) 
are inversely proportional to the boundaries of the power-law 
scaling regime of the duration distribution Ilt(favai))- Results 
agree with MFT predictions (Table |^. 


ulations [39]. The slip events distribute themselves ac¬ 
cording to power laws ~ Ds{S) - 

De{E) ^ PS{u}) ^ with critical exponents 

that are in agreement with the Mean Field interface de- 
pinning transition universality class (Table |l]). 

Our extended results from the PFC model agree with 
the majority of the robust experimental and computa¬ 
tional results in the literature (for an extended sum¬ 
mary see Table |T|. Friedman et al. [21] analyzes the slip 
statistics of compressed crystalline nano-pillars in dif¬ 
ferent stress bins as the flow stress is approached and 
get K = 1.5 (also cr = 2). Dimiduk et al. measure 
K = 1.5 — 1.6 from compression experiments on micro¬ 
pillars at slowly increasing stress [9]. Similarly, in 2D 
DDD and continuum models with quasi-static stress in¬ 
crease in the pinned regime, Zaiser et al. calculate 
AC = 1.4 [HiD]. Slip-event energy amplitudes are power 
law distributed with exponent ke = 1-8 ^ in simula¬ 
tions, and Ke = 1.6 in experiments [1]. Richeton et al. 
[4] report ke = 1-5 for the energy distribution of acous¬ 
tic emission deformation experiments. Of course in [TT] a 
large number of scaling exponents was extracted from 2D 
DDD simulations and as shown in Table [^ they corrobo¬ 
rate that crystalline plasticity belongs to the universality 
class of mean field depinning for the discussed static and 
dynamic properties of the avalanche statistics. This re¬ 
sult is also consistent with analytic results [iniiis] and 
simulations in 3 dimensions m- 
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Discussion: We approach crystalline plasticity with a 
new and alternative simulation methodology: the phase 
field crystal model. The PFC simulations can essen¬ 
tially be thought of as molecular dynamics simulations 
but greatly sped up. At the same time they are free from 
the several phenomenological rules and constraints that 
discrete dislocation dynamics need in order to incorpo¬ 
rate the variety of dynamical phenomena that take place 
in a stressed crystal. The reason is that the phase field 
crystal reproduces faithfully the real crystal including its 
elastic properties and topological behavior. 

By employing this sheared PFC model we were suc¬ 
cessful in extracting the characteristic statistical scaling 
behavior of plastic deformation. We verified the robust¬ 
ness of the DDD simulations, nullified potential artifacts 
of the add-on phenomenological creation and annihila¬ 
tion rules of the DDD simulations and strengthened the 
universal conclusions. Our results reaffirm that all the 
exponents and scaling forms extracted here for crystal 
plasticity are consistent with the mean field interface de- 
pinning nniversality class - even in the absence of frozen- 
in pinning centers. 

An increasing nnmber of studies (including this work) 
have indicated the striking similarities between crystal 
plasticity and the interface depinning dynamic phase 
transition [iniiiiiiis]. The critical exponents we found 
here are in excellent agreement with the mean field the¬ 
ory of the interface depinning universality class (see Table 
|T]) even though the PFC results are performed at finite 
temperature, e ^ ksT, and as a result the extracted scal¬ 
ing relations are plagued by larger fluctuations (due to 
temperature-induced dislocation creep), when compared 
to DDD simulations at T = 0. 
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